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Abstract 


An  unstructured  grid  Large  Eddy  Simulation  (LES)  methodology  hcis  been  developed  for  compressible  high 
speed  flows.  The  filtered  compressible  Navier-Stokes  equations  are  solved  on  an  unstructured  grid  of  tetrahe- 
dra.  The  inviscid  fluxes  are  obtained  from  an  exact  locally  one-dimensional  Riemann  solver  using  Godunov’s 
method.  The  viscous  fluxes  are  obtained  using  a  discrete  analog  of  Gauss’  Theorem,  The  reconstruction 
is  performed  using  a  Least  Squares  technique.  The  temporal  integration  is  a  Runge-Kutta  method.  The 
algorithm  is  overall  second  order  accurate  in  space  and  time.  Four  flowfields  have  been  computed:  decay  of 
isotropic  turbulence,  channel  flow,  supersonic  flat  plate  boundary  layer  and  supersonic  compression  corner. 
The  first  and  second  cases  are  effectively  incompressible,  while  the  third  and  fourth  cases  are  supersonic 
(Mach  3).  The  computed  results  show  close  agreement  with  experiment  and  Direct  Numerical  Simulation, 
and  validate  the  unstructured  grid  LES  methodology. 
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Introduction 


The  effective  design  of  high  speed  aircraft  and  missiles  depends  critically  upon  accurate  prediction  of  aerody¬ 
namic  and  aerothermodynamic  performance  which  are  strongly  affected  by  flow  turbulence  under  most  flight 
conditions.  From  an  engineering  standpoint,  the  aircraft  or  missile  aerodynamicist  needs  the  capability  for 
accurate  prediction  of  the  mean  and  rms  fluctuating  surface  pressure  {pyj  and  and  surface  heat  transfer 
{qyj  and  mean  surface  skin  friction  (f«,),  and  locations  of  primary  and  secondary  separation. 

Table  1  RANS  Capability  for  3-D  Shock  Wave  Boundary  Layer  Interaction 


Quaniiiy 

Saiisfaciory 

Unsaiisfaciory 

No  capabiliiy  shown 

Pw 

V 

Pw 

v 

Qw 

V 

Qw 

yj 

Primary  Separation 

V 

Secondary  Separation 

V 

The  current  methodology  for  prediction  of  compressible  turbulent  flows  is  based  on  the  Reynolds-averaged 
Navier-Stokes  (RANS)  equations  (Knight  1993).  This  approach  has  yielded  a  hierarchy  of  turbulence  models 
extending,  from  zero-equation  to  full  Reynolds  Stress  Equation  models.  While  these  models  have  generally 
been  capable  of  predicting  the  engineering  quantities  of  interest  in  weakly  perturbed  boundary  layers,  they 
have  been  unable  to  accurately  predict  the  complex  3-D  flows  which  are  encountered  in  highly  maneuvering, 
high  angle-of-attack  flight.  Two  recent  extensive  reviews  have  documented  the  capabilities  and  deficiencies 
of  a  wide  range  of  RANS  models  for  prediction  of  complex  3-D  flows  with  shock  wave-turbulent  boundary 
layer  interactions  (Knight  1997,  Knight  and  Degrez  1998).  The  results,  summarized  in  Table  1,  indicate 
that  a  significant  number  of  critical  engineering  quantities  are  not  capable  of  prediction  by  current  RANS 
models.  Therefore,  more  advanced  turbulence  models  are  needed  which  have  the  ability  to  simulate  the 
complex  physics  of  turbulence  with  greater  generality. 

Large  Eddy  Simulation  (LES)  is  an  alternative  to  RANS  hich  may  be  capable  of  predicting  more  (or  all) 
of  the  aerodynamic  and  aerothermodynamic  quantities  of  engineering  interest  described  above.  In  LES,  the 
governing  equations  are  spatially  filtered  on  the  scale  of  the  numerical  grid.  The  large,  energy-containing 
eddies  are  directly  computed.  These  eddies  are  strongly  influenced  by  the  physical  geometry  and  configura¬ 
tion  of  the  flow.  Thus,  the  direct  computation  of  the  large  eddies  by  LES,  as  opposed  to  the  modeling  of  the 
large  eddies  by  RANS,  gives  greater  generality,  in  principle,  to  LES.  The  influence  of  the  unresolved  scales 
of  motion  is  simulated  using  a  subgrid-scale  (SGS)  model  (Smagorinsky  1963,  Lilly  1967,  Deardorff  1970, 
Germane  ei  al  1991,  Piomelli  et  al  1991,  Ghosal  ei  al  1995)  or  by  the  inherent  dissipation  in  the  numerical 
scheme  (Boris  ei  al  1992,  Oran  and  Boris  1993,  Porter  ei  al  1994,  Grinstein  1996,  Ansari  and  Strang  1996). 
Because  the  statistics  of  the  small  scale  turbulence  are  expected  to  be  more  homogeneous  and  isotropic  than 
those  of  the  large  scales,  a  general  model  of  the  small  scales  seems  more  plausible  than  a  general  model  of 
the  entire  spectrum  of  turbulent  motions. 

LES  has  been  shown  to  be  both  a  useful  research  tool  for  understanding  the  physics  of  turbulence,  and 
also  a  predictive  method  for  flows  of  engineering  interest.  Recent  compendia  and  reviews  include  Galperin 
and  Orszag  (1993),  Mason  (1994),  Lesieur  and  Metais  (1996)  and  Moin  (1997).  Many  models  have  been 
developed  for  the  subgrid-scale  stress  tensor.  These  include  the  conventional  Smagorinsky  eddy  viscosity 
model  (Smagorinsky  1963,  Lilly  1967,  Deardorff  1970),  the  spectra  eddy  viscosity  model  of  Kraichnan  (1976), 
the  dynamic  SGS  model  of  Germane  ei  al  (1991),  the  scale  similarity  model  of  Bardina  ei  al  (1980),  and  the 
localized  dynamic  SGS  model  of  Ghosal  ei  al  (1995)  and  more  recently  of  Menon  and  Kim  (1996),  and  many 
others.  Although  most  research  has  focused  on  incompressible  turbulent  flows,  there  has  recently  emerged  a 
growing  interest  in  applications  of  LES  to  compressible  turbulent  flows.  Examples  include  Yoshizawa  (1986), 
Speziale  ei  al  (1988),  Moin  ei  al  (1991),  Erlebacher  ei  al  (1992),  Zang  ei  al  (1992),  El-Hady  ei  al  (1994), 
Jansen  (1997),  Spyropoulos  and  Blaisdell  (1996),  and  Haworth  and  Jansen  (1996).  Nearly  all  compressible 
LES  has  employed  spectral  methods  or  structured  grids,  with  the  exception  of  Jansen  and  Haworth. 
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Apart  from  the  complexities  of  the  flowfield,  the  complicated  geometries  of  high  speed  vehicles  is  also  a 
challenge.  To  enable  treatment  of  complex  geometries  and  also  achieve  high  resolution  of  the  flowfield 
dynamically,  we  employ  an  unstructured  grid.  There  are  two  important  advantages  of  unstructured  grids. 
First,  algorithms  have  been  developed  to  facilitate  automatic  generation  of  unstructured  grids  for  a  complex 
geometries  (see,  for  example,  the  discussion  in  Barth  (1990,  1992).  These  grid  generation  methods  can  be 
substantially  more  efficient  (in  terms  of  user  time)  than  some  of  the  multi-block  structured  grid  generation 
methods  used.  Second,  local  mesh  refinement,  either  adaptive  or  fixed,  can  been  performed  much  more 
readily  for  unstructured  grids. 

The  report  summarizes  the  research  in  Large  Eddy  Simulation  of  compressible  turbulent  flows  using  unstruc¬ 
tured  grids.  Two  methods  for  simulation  of  the  subgrid  scale  stresses  have  been  examined.  The  first  method 
is  the  Monotone  Integrated  Large  Eddy  Simulation  (MILES)  technique.  The  second  method  is  a  hybrid 
technique  combining  MILES  with  a  Smagorinsky  eddy  viscosity  model  for  the  subgrid  scale  stresses.  These 
two  methods,  together  with  different  algorithms  for  the  inviscid  fluxes  and  function  reconstruction,  have  been 
evaluated  for  four  turbulent  flows:  the  decay  of  isotropic  incompressible  turbulence,  channel  flow,  supersonic 
boundary  layer  and  supersonic  compression  corner.  The  results  are  in  overall  good  agreement  with  the 
experiment  and  Direct  Numerical  Simulation  (DNS),  thereby  validating  the  accuracy  of  the  methodology. 

Governing  Equations 


The  governing  equations  are  the  three-dimensional  filtered  Navier-Stokes  equations.  For  a  function  /,  its 
filtered  form  /  is 


/ 


GfdV 


where  G  is  the  filtering  function,  and  its  Favre-averaged  form  /  is 

P 

where  p  is  the  density.  From  the  Navier-Stokes  equations  for  the  instantaneous  flow  variables  density  (/>), 
velocity  in  the  ith  coordinate  direction  (w,),  pressure  (p)  and  temperature  (T),  Favre-averaging  and  spatial 
filtering  yield  the  filtered  Navier-Stokes  equations  (here  written  using  the  Einstein  summation  notation  where 
repeated  indices  denote  summation) 


dt  dxi,  ~ 
dpUi  ^  dpUiUk 


dpi 


dt 

d 


dxk 


dp  dTik 

dxi  dxk 


dOu  d 

+  _(^a-  +  fl«  =  ^  +  — (3;«) 


p  =  pRT 


where 


%k  =  Tik  +  (Tik 

Tik  =  -p{uiUk  ■ 


UiUk) 

duj 


2duj  dui  duk 

3'^^ik  +  - 1- 


dik  =  pif)  (- 

Qk  =  Qk+  Qk 

Qk  =  —pcp{Tuk  —  TUk) 

~  df 

pi  =  pc„T  +  +  pk 

pk  =  §(ptliu7~ 


dxk  dx 


7) 
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Two  different  Sub-Grid-Scale  (SGS)  models  are  employed.  The  first  model  is  Monotone  Integrated  Large 
Eddy  Simulation  (MILES)  wherein  the  numerical  algorithm  itself  provides  the  requisite  dissipation  associated 
with  the  subgrid  scale  motions.  The  second  model  is  the  classical  constant-coefficient  Smagorinsky  method 

"  2  [q^.  +  QiJ 
'^ij  —  2Cjip^  ^ Smn  Smn  (^^ij 

Qj  =  ^ 

where  Cr  =  0.00423  and  A  is  the  length  scale  which  is  related  to  the  local  grid  size.  For  boundary  layer 
flows,  A  is  multiplied  by  the  Van  Driest  damping  factor 

D  =  l- 

where  A  =  26,  =  nurjvw  is  the  normal  distance  to  the  (nearest)  solid  boundary  normalized  by  the 

viscous  length  scale  Vy^jur  where  is  the  kinematic  viscosity  evaluated  at  the  wall  and  XJr  is  the  local 
friction  velocity. 

We  simplify  the  notation  by  hereafter  dropping  the  tilde  ^  and  overbar  “.  The  flow  variables  are  nondi- 
mensionalized  using  the  reference  density  velocity  {/oo?  static  temperature  Too  and  length  scale  L,  with 
Mach  number  Mq©  =  Uool\/lRToo>  The  governing  equations  are  therefore 

dp  dpuk  _  - 
dt  dxi: 

dpuj  dpUjUj:  __  ,  dTik 

dt  dxk  dxi  dxk 


P  = 


9xk 

pT 

'iMl 


Numerical  Algorithm 

The  governing  equations  are  expressed  in  finite  volume  form  for  a  control  volume  V  with  surface  dV 

4  /  QdV+  f  iFt  +  Gj+Hk)-ndA  =  0 

at  Jy  Jqy 

where  Q  is  the  vector  of  dependent  variables 


and  the  flux  vectors  are 


F=  ( 


'  pu 

'  ) 

'  pv) 

pu^  +p-  Tix 

pUV  -  Txy 

puw  -  Txz 

PUV  -  %y 

G=< 

pV^  +P-Tyy  >  ,  H  =  i 

pVW  -  Tyz 

puw  —  Txz 

pVW  -  Tyz 

pw^  +  p-  Tzz 

.  (pe  +  p)u  -  Qx-Px  ^ 

.  {pe  +  P)V  -Qy-^y  ) 

.  {pe  +  p)w  -Qz-Pz  ^ 
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with 


Px  4"  4" 

—  ^xy'^  4*  ^yy*^  4”  ^yz"^ 

Pz  —  '^xz'^  4-  %fz'^  4-  TzzW 

An  unstructured  grid  of  tetrahedra  is  employed,  with  a  cell-centered  storage  architecture.  The  cell-averaged 
values,  stored  at  the  centroid  of  each  tetrahedron  of  volume  Vi  are 

=  ^  !  Qdv 

M  JVi 

The  inviscid  fluxes  are  computed  using  Godunov’s  method  which  is  an  exact  one-dimensional  Riemann  solver 
(Gottlieb  and  Groth  1988)  applied  normal  to  each  face.  The  inviscid  flux  computations  require  the  values 
of  each  variable  on  either  side  of  the  cell  faces.  These  values  are  obtained  from  the  cell-averaged  values 
by  second-order  or  third-order  function  reconstruction  using  the  Least  Squares  method  of  Ollivier-Gooch 
(Ollivier-Gooch  1997).  The  second-order  function  reconstruction  method  of  Prink  (1994)  was  employed  in 
some  of  the  earlier  LES  studies,  but  was  found  inferior  to  the  method  of  Ollivier-Gooch  (Okong’o  and  Knight 
1998).  More  details  on  the  reconstruction  schemes  are  given  in  Okong’o  and  Knight  (1998). 

The  viscous  fluxes  and  heat  transfer  are  computed  by  application  of  Gauss’  theorem  to  the  control  volume 
whose  vertices  are  the  centroids  of  the  cells  which  share  each  node.  The  second-order  accurate  scheme  (in 
2-D)  is  given  by  Knight  (1994)  and  the  extension  to  3-D  is  straightforward. 

Parallelization 

The  code  is  parallelized  using  domain  decomposition  and 
Message  Passing  Interface  (MPI).  Domain  decomposition  is 
performed  in  a  pre-processing  step.  The  domain  is  decom¬ 
posed  in  a  single  direction  with  equal  number  of  tetrahedra 
in  each  domain.  A  halo  of  cells  is  added  in  each  domain 
to  provide  data  on  the  adjacent  domain,  and  the  halo  cell 
data  is  updated  at  every  subiterate  of  the  time  integration. 

An  example  is  shown  in  Fig.  1  for  the  LES  of  decay  of 
isotropic  turbulence.  The  numerical  algorithm  achieves  ex¬ 
cellent  parallel  performance.  For  example,  the  speed-up  on 
four  processors  of  the  SGI  Power  Onyx  with  R- 10000  pro- 
cessors  is  3.7  for  93%  efficiency  (Knight  et  al  1998).  pig  j  Example  of  domain  decomposition 

Results 

Four  different  configurations  have  been  examined:  decay  of  isotropic  turbulence,  channel  flow,  flat  plate 
boundary  layer  and  supersonic  compression  corner.  The  first  and  second  cases  are  effectively  incompressible, 
while  the  third  and  fourth  cases  are  supersonic  (Mach  3). 

Decay  of  Isotropic  Turbulence 

The  benchmark  experiment  by  Comte-Bellot  and  Corrsin  (CBC)  (1971)  of  the  decay  of  isotropic  turbulence 
was  selected  to  examine  the  following  specific  issues  of  the  unstructured  grid  LES  code:  the  MILES  vs 
hybrid  SGS  models,  the  method  for  reconstruction  of  the  flow  variables  to  the  cell  faces,  the  demonstration 
of  achieving  consistent  results  using  two  different  grids,  the  type  of  Riemann  solver,  regular  vs  random  grids, 
the  convergence  parameter  in  Godunov’s  method,  and  the  number  of  Fourier  intervals  used  in  the  initial 
spectrum.  In  the  experiment,  turbulence  was  generated  using  a  biplane,  square  rod  grid  with  mesh  size 
M  =  5.08  cm  and  solidity  of  0.34  in  a  uniform  mean  flow  of  velocity  Uq  =  10  m/s.  The  Reynolds  number 
based  on  the  grid  spacing  is  Re  =  34000.  The  measurements  were  performed  downstream  of  the  grid  at  three 
locations,  Uot/M  =  42,  98,  and  171  where  t  is  the  dimensional  ‘‘time”  in  the  experiment.  The  computation 
was  performed  in  a  periodic  box.  Details  are  presented  in  Knight  et  al  (1998). 
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Ten  different  computations,  as  summarized  in  Table  2,  were  performed  to  evaluate  the  issues  described 
below.  Cases  1  through  7  and  10  use  the  MILES  model  (Le.,  no  SGS  model),  while  Ccises  8  and  9  use  the 
hybrid  model  (i.e.,  MILES  plus  the  Smagorinsky  SGS  model).  Case  1  employs  a)  second-order  least-squares 
reconstruction  with  8  cells  in  the  stencil,  b)  a  regular  grid,  c)  Godunov’s  method,  d)  a  convergence  parameter 
of  10”^^  in  Godunov’s  method,  and  e)  Nk  =  26  in  the  initial  energy  spectrum.  Cases  2  through  7  in  succession 
modify  one  of  these  five  choices.  For  example.  Case  2  employs  Frink’s  method  for  reconstruction.  Cases  8 
and  9  use  the  hybrid  (MILES  plus  Smagorinsky  models)  with  different  length  scale  A.  The  Smagorinsky 
constant  Cr  is  chosen  to  be  0.012.  Cases  1  to  9  employ  Grid  1,  and  Case  10  employs  Grid  2. 


Table  2  Decay  of  Isotropic  Turbulence 


Case 

Gnd 

SGS 

A 

Reconstruction 

Grid 

Ni 

Riemann 

Toler 

Model 

Method 

Type 

Solver 

1 

1 

MILES 

n/a 

LS  2nd 

Reg 

26 

Godunov 

10-10 

2 

1 

MILES 

n/a 

F  2nd 

Reg 

26 

Godunov 

10-10 

3 

1 

MILES 

n/a 

F2nd 

Reg 

26 

Roe 

4 

1 

MILES 

n/a 

First 

Reg 

26 

Godunov 

10-10 

5 

1 

MILES 

n/a 

F  2nd 

Ran 

26 

Godunov 

10-10 

6 

1 

MILES 

n/a 

F  2nd 

Reg 

26 

Godunov 

10-12 

7 

1 

MILES 

n/a 

F  2nd 

Reg 

52 

Godunov 

10-10 

8 

1 

Hybrid 

F  2nd 

Reg 

26 

Godunov 

10-10 

9 

1 

Hybrid 

F  2nd 

Reg 

26 

Godunov 

10-10 

10 

2 

MILES 

n/a 

F  2nd 

Reg 

26 

Godunov 

10-10 

Legend 


Grid  1  No.  of  tetrahedra  =  163,840 

Grid  2  No.  of  tetrahedra  =  1,310,720 

MILES  Monotone  Integrated  Large  Eddy  Simulation  (Boris  et  al  1992) 

LS  2nd  Second  order  least-squares  reconstruction  method 

F  2nd  Second  order  Frink’s  reconstruction  method 

First  First  order  reconstruction  method 

Reg  Regular  grid  of  tetrahedra 

Ran  Grid  obtained  by  random  perturbing  the  nodes  of  the  regular  grid 

Nk  Number  of  Fourier  intervals  used  for  initial  energy  spectrum 

Toler  Tolerance  employed  in  iteration  solution  for  p* 


The  decay  of  the  resolved  turbulence  kinetic 
energy  for  Cases  1  to  9  is  compared  with  the 
filtered  experimental  data  of  CBC  in  Fig.  2. 

The  MILES  method  using  either  second  order 
reconstruction  method  (Cases  1  to  3,  and  5  to  | 

7)  accurately  predicts  the  turbulence  decay.  ^ 

This  remarkable  result  implies  that  the  inher-  | 

ent  numerical  dissipation  in  the  numerical  al-  8 

gorithm  (due  to  the  finite  order  accuracy  of  | 

the  reconstruction  and  quadrature)  provides 
a  reasonable  model  of  turbulent  energy  dissi¬ 
pation.  The  computations  employing  the  hy¬ 
brid  (MILES  plus  Smagorinsky  models)  ex¬ 
hibit  only  a  small  difference  compared  to  the 
MILES  simulations,  since  the  dissipation  of 
turbulence  energy  is  almost  entirely  a  conse¬ 
quence  of  the  numerical  algorithm  (MILES). 


Fig.  2  Filtered  turbulent  kinetic  energy 
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The  computation  using  the  second  order  least-squares  method  (Case  1)  shows  a  lower  decay  rate  than 
those  using  Frink’s  method  (Cases  2,  3,  and  5  to  7).  The  second  order  least-squares  evidenlty  has  a  lower 
inherent  numerical  dissipation  than  the  second  order  Frink’s  reconstruction.  An  SGS  model  may  be  needed 
to  compensate  the  reduced  inherent  numerical  dissipation  in  order  to  achieve  the  best  agreement  with 
experiment.  The  computation  using  the  first  order  reconstruction  (Ceise  4)  shows  poor  agreement  with 
experiment  as  anticipated. 

The  MILES  approach  yields  consistent  results  when  the  grid  is  refined  (Knight  ei  al  1998).  The  results  using 
either  Roe’s  method  (Case  3)  or  Godunov’s  method  (Case  2)  are  essentially  identical.  This  is  expected  since 
the  velocity  fluctuations  are  small  compared  to  the  speed  of  sound. 

The  choice  of  regular  vs  random  grids  ha^  no  significant  effect  on  the  accuracy  of  the  computation  as  seen  in 
Fig.  2.  where  the  turbulence  energy  for  Cases  2  and  5  are  essentially  identical.  Moreover,  the  variation  in 
the  convergence  parameter  used  in  Godunov’s  method  (Cases  2  and  6),  and  the  number  of  Fourier  intervals 
in  the  initial  condition  (Cases  2  and  7)  have  negligible  effect. 

a)  b) 


Fig.  3  Energy  spectrum  at  Uot/M  =  98  Fig.  4  Energy  spectrum  at  UotfM  =  171 

The  effect  of  the  functional  reconstruction  method  on  the  turbulence  energy  spectrum  is  shown  in  Figs.  3  and 
4  where  the  experimental  data  and  computed  results  using  the  second  order  least-squares  and  Frink’s  method 
are  shown  at  UotfM  =  98  and  171,  respectively,  for  Grid  1.  The  least-squares  reconstruction  is  observed  to 
be  more  accurate  than  Frink’s  method.  Since  the  CPU  time  for  the  two  methods  is  approximately  the  same, 
the  least-squares  reconstruction  is  preferable.  A  ‘‘pile-up”  of  energy  at  lower  wave  numbers  is  noted.  This 
wcLS  also  observed  by  Haworth  and  Jansen  (1996)  who  employed  a  similar  grid  resolution. 

Channel  Flow 

The  Reynolds  number  based  on  the  channel  height  and  bulk  velocity  is  5600.  The  Mach  number  is  0.5; 
however,  the  static  temperature  variation  across  the  channel  is  small  (less  than  4%)  so  the  flow  is  effectively 
incompressible.  The  mean  flow  is  in  the  x-  direction,  with  y  and  z  representing  wall  normal  and  span  wise 
directions,  respectively. 

Computations  have  been  performed  using  the  constant-coefficient  Smagorinsky  sub-grid  scale  model  (Cr  = 
0.012,  Prt=  0.4).  The  coefficient  Cr  is  multiplied  by  the  van  Driest  damping  factor  1  —  integrate 

to  the  walls.  The  length  scale  A  for  the  SGS  model  is  the  nodal  spacing  in  the  y-direction.  A  body-force 
is  applied  in  this  direction  and  is  adjusted  at  each  timestep  to  keep  the  bulk  velocity  constant.  Using  the 
global  x-momentum  conservation  for  fully  developed  flow,  fi  =  — 2uJ.  The  bulk  velocity  is  also  the  reference 
velocity  and  so  is  set  equal  to  1.  The  boundary  conditions  are  periodic  in  the  streamwise  (x)  and  spanwise 
(z)  directions  with  no-slip  boundaries  in  the  wall-normal  (y)  direction,  with  isothermal  walls. 

The  grid  size  is  27r  x  1  x  2x/3,  where  the  lengths  have  been  non-dimensionalized  by  the  channel  height.  The 
grid  has  65  x  65  x  65  nodes  (274,625  nodes),  1,310,720  cells,  and  2,646,016  faces.  The  grid  is  stretched  in  the 
y-direction,  with  minimum  spacing  (at  the  walls)  of  0.00278  and  a  maximum  spacing  (in  the  channel  center) 
of  0,0457.  The  initial  condition  is  obtained  from  interpolating  the  fully-developed  solution  on  a  grid  with 
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33  X  65  X  33  nodes  with  the  same  wall-normal  (t/)  spacing  but  double  the  streamwise  and  spanwise  spacing. 
The  flow  has  been  simulated  for  about  twenty-two  flow-through  times,  with  statistics  accumulated  over  the 
last  twelve  flow-through  times.  (One  flow-through  time  is  the  domain  length  in  the  ai-direction  divided  by 
the  bulk  velocity,  i.e.,  27r.)  The  results  are  compared  with  the  experimental  results  of  Eckelmann  (1974)  and 
Kreplin  and  Eckelmann  (1974)  and  Direct  Numerical  Simulation  of  Kim  ei  al  (1987).  In  the  presentation  of 
results,  statistics  are  accumulated  at  all  grid  points.  Planar  averages  are  obtained  by  averaging  at  all  points 
on  planes  parallel  to  the  walls,  thus  making  the  statistics  functions  of  the  wall-normal  (t/)  direction  only. 

The  computed  friction  velocity,  normalized  by  the  bulk  velocity,  is  Ur  =  0.0594  which  is  within  8%  of  the 
DNS  and  experimental  values  of  0.0643.  The  planar  average  of  the  time-averaged  streamwise  velocity  U  is 
shown  in  Figs.  5  and  6.  The  centerline  velocity  is  1.155  which  is  within  0.6%  of  the  experimental  value  of 
1.162.  The  LES  is  close  to  the  DNS  in  the  near- wall  region,  but  overshoots  it  slightly  in  the  core  of  the 
channel,  consistent  with  the  lower  shear  stress  of  the  LES.  Notably,  the  LES  h^ls  excellent  agreement  with 
the  experimental  data  of  Eckelmann  (1979),  whose  log-law  of  =  2.65  In -f  5.9  is  based  on  data  at  Rcc 
of  5600  and  8200.  The  LES  Rcc  of  about  6470  is  between  these  values.  The  experimental  data  points  shown 
in  Fig.  6  are  for  Rcc  —  5600. 


Y  y* 

Fig.  5  Mean  velocity  in  entire  channel  Fig.  6  Mean  velocity  in  wall  units 


Figs.  7,  8  and  9  compare  the  root-mean-square  velocity  fluctuations  from  the  LES  with  the  DNS  results  of 
Kim  et  al  and  the  experimental  results  from  Kreplin  and  Eckelmann  (1979).  In  the  experiment,  Rcc  =  7700 
and  Rcr  =  388,  which  are  slightly  higher  that  the  DNS  Rcc  =  6600  and  Rcr  =  360  and  the  LES  Rcc  ^ 
6470  and  Rcr  =  330.  The  velocity  fluctuations  from  the  LES,  DNS  and  experiment  show  generally  good 
agreement.  The  streamwise  fluctuations  (Fig.  7)  of  the  LES  compare  favorably  with  the  experiment,  but  are 
slightly  higher  than  the  DNS.  The  LES  wall-normal  fluctuations  to  be  lower  than  those  of  the  experiment  and 
DNS  (Fig.  8);  however,  the  DNS  is  also  significantly  lower  than  the  experiment.  The  spanwise  fluctuations 
(Fig.  9)  show  comparable  agreement  with  experiment  for  the  LES  and  DNS.  The  computed  Reynolds  stress, 
shown  in  Fig.  10,  compares  very  well  with  the  DNS. 


Fig.  7  Rms  streamwise  velocity  fluctuations  Fig.  8  Rms  normal  velocity  fluctuations 
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Fig.  9  Rms  spanwise  velocity  fluctuations  Fig.  10  Reynolds  stress 

The  higher-order  velocity  statistics  based  on  the  resolved  velocities  are  shown  in  Figs.  11  to  16.  The  skewness 
for  the  streamwise  velocity  is  defined  by  5(u')  =  u*u*u* and  likewise  for  and  ly',  where  the  bar 
denotes  the  time-average  and  the  prime  denotes  fluctuation  relative  to  the  time-average.  Fig.  11  shows  the 
streamwise  skewness  from  the  LES  to  be  larger  than  the  DNS,  particularly  in  the  center  of  the  channel.  The 
reaisons  for  the  disagreement  with  DNS  are  under  evaluation.  However,  the  agreement  for  the  wall-normcil 
and  spanwise  skewness,  Figs.  12  and  13,  is  excellent.  The  flatness  for  the  streamwise  velocity  is  defined  by 
F(u')  =  u'w'ii'tz/w'ti'  and  likewise  for  v'  and  w\  The  streamwise  statistics  of  the  LES  deviate  from  the 
DNS  (Fig.  14),  while  the  normal  and  spanwise  flatness  profiles,  shown  in  Figs.  15  and  16,  agree  very  well 
with  the  DNS. 


Fig.  11  Streamwise  velocity  fluctuation  skewness  Fig.  12  Normal  velocity  fluctuation  skewness 


Fig.  13  Spanwise  velocity  fluctuation  skewness  Fig.  14  Streamwise  velocity  fluctuation  flatness 
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Fig.  15  Normal  velocity  fluctuation  flatness 


Fig.  16  Spanwise  velocity  fluctuation  flatness 


The  channel  flow  LES  simulation  is  effectively  MILES,  im¬ 
plying  that  the  SGS  dissipation  is  attributable  to  the  nu¬ 
merical  algorithm.  This  is  evident  in  Fig.  17  which  shows 
the  Smagorinsky  eddy  viscosity  normalized  by  the  molecu¬ 
lar  viscosity.  The  peak  value  of  the  Smagorinsky  is  less  than 
10%  of  the  molecular  viscosity  and  occurs  in  the  middle  of 
the  channel.  Thus,  the  Smagorinsky  eddy  viscosity  does 
not  significantly  contribute  to  the  energy  dissipation  which 
is  due  almost  entirely  to  the  numerical  algorithm.  This  is 
a  particularly  important  result  and  the  first  observation  of 
this  type  for  an  unstructured  grid  computation  of  turbulent 
channel  flow. 


It  is  recognized  that  the  above  comparison  between  the  DNS  and  the  LES  is  not  precise  because  the  DNS 
has  not  been  filtered  onto  the  LES  grid.  However,  for  this  case  most  of  the  energy  is  in  the  resolved  scales, 
as  shown  by  Fig.  17. 

Supersonic  Flat  Plate  Turbulent  Boundary  Layer 

An  adiabatic  flat  plate  turbulent  boundary  layer  at  Mach  3  and  Reynolds  number  Res  =  2  x  10^  (based 
on  the  incoming  boundary  layer  thickness  S)  has  been  computed.  The  Reynolds  number  based  on  the 
momentum  thickness  S2  and  wall  viscosity  is  Res^  =  600.  The  Reynolds  number  is  sufficently  high  to 
achieve  turbulent  flow. 


Fig.  18  Computational  domain  Fig.  19  Instaneous  plot  of  T 

The  inflow  conditions  are  obtained  using  a  compressible  extension  of  the  method  of  Lund  et  al  (1998).  The 
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simulation  generates  its  own  inflow  conditions  through  a  sequence  of  operations  where  the  velocity  field  at 
a  downstream  station  is  rescaled  and  reintroduced  at  the  inflow  boundary  (Fig.  18).  Defining  x^y  and  z 
to  denote  the  streamwise,  transverse  and  spanwise  directions,  respectively,  the  size  of  the  computational 
domain  is  Lx  =  14. 8^,  Ly  =  3.4<5  and  Lz  =  4.4^.  The  spanwise  width  Lz  is  approximately  six  times  the 
experimental  spanwise  streak  spacing  (assuming  the  compressible  turbulent  boundary  layer  streaks  scale 
in  accordance  with  incompressible  experimental  results).  The  streamwise  length  Lx  is  approximately  three 
times  the  mean  experimental  streamwise  streak  size.  The  height  Ly  is  based  on  the  requirement  that  acoustic 
disturbances  originating  at  the  upper  boundary  do  not  interact  with  the  boundary  layer  on  the  lower  wall. 
The  grid  resolution  is  1,600,000  tetrahedral  cells.  The  grid  at  the  wall  uses  101  x  65  nodes.  The  nodal  grid 
is  stretched  uniformly  in  the  y  using  a  geometric  factor  of  1.088.  The  height  An  of  the  first  cell  adjacent 
to  the  boundary  is  less  than  one  wall  unit  (Anur/uw  <  1  where  i/y,  is  the  kinematic  viscosity  at  the  wall, 
=  y/^wIPw  is  the  friction  velocity,  Ty,  is  the  wall  shear  stress  and  is  the  density  at  the  wall).  The 
initial  condition  is  a  turbulent  mean  profile  with  random  fluctuations.  The  simulation  is  run  first  for  90 
inertial  timescales  6/Uoo  in  order  to  eliminate  starting  transients  (Lund  ei  al  1998). 

Fig.  19  shows  the  instantaneous  temperature  contours  within  the  boundary  layer.  Large  turbulent  streak-like 
structures  eject  high  temperature  fluid  outward  the  inner  layer. 

For  a  function  /,  its  average  in  time  form  <  /  >  is  defined  by 

<f>=; — r/ 

V  “  Jti 


and  its  time  fluctuating  part  is 


/"  =  /-</> 


In  order  to  provide  converged  data,  the  primitive  variables  are  averaged  in  spanwise  direction  and  the 
statistical  evaluations  are  perfomed  on  a  period  longer  than  tf  —  ti  =  406/Cfcx>*  The  notation  for  the 
combined  temporal  and  spanwise  average  is 


<f>= 


1  1  r^‘ 

~  Jo  Jti 


fdtdz 


A  simplifying  notation  is  used  for  the  velocity,  temperature  and  pressure 


U  =<^u> 


The  predicted  friction  velocity  is  Ur/Uoo  —  0.053  (the  friction  coefficient  is  Cf  =  2.10  x  10“^)  which  is  within 
3%  of  the  theoretical  value  of  0.0544  b2Lsed  on  the  friction  law  obtained  from  the  combined  Law  of  the  Wall 
and  Wake  (using  the  classical  coefficients  /c  =  0.41;  (7  =  5.0;  11  =  0.55).  The  uncertainty  in  the  experimental 
data  of  skin  friction  is  typically  ±5%,  and  the  accuracy  of  the  correlation  is  typically  ±10%  as  discussed  by 
Hopkins  and  Inouye  (1971). 

The  computed  adiabatic  wall  temperature  is  Taw/Too  =  2.68  which  is  within  3%  of  the  theoretical  value  of 
'Law /Too  =  2.602  obtained  from  the  empirical  formula 


Ml 


where  Prt^  =  0.89  is  the  mean  turbulent  Prandtl  number. 
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Fig.  20  Mean  temperature  vs  velocity  Fig.  21  Mean  velocity 

In  Fig.  20  the  temperature-velocity  relationship  is  evaluated.  The  data  fit  closely  to  the  classical  quadratic 
profile  originally  due  to  Walz  (1969). 

The  log-linear  plot  of  the  mean  streamwise  velocity  profile  is  shown  in  Fig.  21.  The  velocity  is  rescaled 
according  to  the  Van-Driest  transformation  and  the  wall  friction  velocity.  The  height  y'^  is  in  wall  units. 
The  logarithmic  regions  appears  in  the  range  30  <  t/*^  <  150.  In  this  representation  the  velocity  appears  to  be 
over-predicted  when  compared  with  the  empirical  and  classical  approximation  UvdIUt  =  2.51og(i/'^)  -h  5.1. 
The  predicted  profile  fits  more  accurately  with  the  equation  UvulUr  =  2.91og(2/‘^)  +  5.7  (dashed  line).  In 
fact,  the  difference  is  not  as  large  than  it  first  appears.  The  representation  is  very  sensitive  with  the  exact 
friction  velocity  value  {y  is  measured  in  wall  units,  which  is  wall  friction  function).  If  Ur  is  assumed  to  be 
equal  to  0.0598  then  the  logarithm  region  will  fit  perfectly  with  the  classical  equation.  One  explanation  of 
this  difference  could  be  the  inflow  regeneration  method  that  needs  to  be  adjusted  accurately  (Lund  ei  al 
1998).  Further  tests  are  in  progress  in  order  to  evaluate  the  convergence  of  the  constants  2.9  and  5.7  with 
mesh  refinement  as  well. 

Table  3  Flat  Plate  Boundary  Layer  Experimental  Data 


Name 

Mach  No. 

Res 

LES 

3.0 

20  X  10® 

DNS  Adams  (1997) 

3.0 

25  X  10® 

Johnson  &  Rose  (1975) 

2.9 

1000  X  10® 

Konrad  (1993) 

2.9 

1590  X  10® 

Konrad  &  Smits  (1998) 

2.87 

1900  X  10® 

Muck  ei  al  (1984,  1985) 

2.87 

1638  X  10® 

Zheltovodov  ei  al  (1986) 

1.7-9.4 

up  to  2000  X  10® 

All  cases  eire  adiabatic  wall 


The  mean  streamwise  resolved  turbulent  kinematic  normal  stress  normalized  using  the  local  mean 

density  <C  p  >  and  wall  shear  stress  is  shown  in  Fig,  22.  As  discussed  in  Zheltovodov  and  Yakovlev 
(1986)  and  Smits  and  Dussauge  (1996),  the  scaling  <C  p  /r^;  provides  an  approximate  self¬ 

similar  correlation  of  experimental  data  for  supersonic  flat  plate  zero  pressure  gradient  adiabatic  boundary 
layers,  although  the  meeisurements  close  to  the  wall  are  subject  to  considerable  uncertainty.  In  Fig.  22 
data  are  displayed  from  Konrad  and  Smits  (1998),  Johnson  and  Rose  (1975),  Muck  ei  al  (1984,  1985), 
Konrad  (1993),  as  well  as  upper  and  lower  bounds  of  an  extensive  set  of  experimental  data  for  the  Mach 

number  range  M  —  1.72 - 9.4  in  accordance  with  generalizations  of  Zheltovodov  and  Yakovlev  (1986). 

The  characteristics  of  the  different  experiments  are  displayed  in  Table  3.  The  computed  results  show  good 
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agreement  with  experiment  for  the  main  part  of  the  boundary  layer  {y/8  >  0.2),  despite  a  significantly 
higher  experimental  Reynolds  number.  The  decreasing  slope  corresponds  precisely  to  Johnson  and  Rose 
(1975)  data.  For  yfS  <  0.2  the  presence  of  the  typical  high  level  peak  in  the  near  wall  region  is  supported  by 
experimental  data  of  Konrad  (1993)  and  the  Direct  Numerical  Simulation  data  from  Adams  (1997),  which 
is  nearly  at  the  same  Reynolds  number  as  the  LES.  However,  no  conclusion  can  be  drawn  about  the  precise 
y  position  and  the  width  of  this  peak  without  further  experimental  data  or  DNS. 


Fig.  22  Reynolds  normal  stress  Fig.  23  Reynolds  shear  stress 

In  Fig.  23  the  Reynolds  shear  stress  distributions  are  shown  for  the  same  experiments  and  the  DNS.  Again, 
the  data  fit  well  in  the  outer  part  of  the  boundary  layer.  The  maximum  value  and  the  decreasing  slope  are 
again  well  predicted. 

Compression  Corner 

Supersonic  flow  past  a  compression  corner  is  an  important  problem  in  aerodynamics.  It  represents,  for 
example,  the  deflection  of  a  control  surface  on  a  wing.  The  shock  can  cause  a  boundary  layer  separation 
upstream  of  the  point  of  impingement  of  the  primary  shock,  with  a  secondary  shock  forming  near  the 
separation  (Andreopoulos  and  Muck  1987,  Dolling  and  Or  1983,  Horstman  ei  al  1977,  Settles  et  al  1979,  Smits 
and  Muck  1987,  Zheltovodov  et  al  1983,  Zheltovodov  and  Yakolev  1986,  and  Zheltovodov  1996)  Reynolds- 
averaged  Navier-Stokes  simulations  have  failed  to  accurately  predict  the  flow  characteristics  (Knight  and 
Degrez  1998)  such  as  fluctuating  pressure  and  heat  transfer. 

Table  4  Compression  Corner  Experimental  Data 


Name 

a 

Mach 

Res 

Res2 

LES 

8" 

3.0 

20 

X 

10® 

600 

DNS  Adams  (1997) 

18° 

3.0 

25 

X 

10^ 

n/a 

Muck  et  al  (1984) 

8° 

2.87 

1638 

X 

10^ 

40800 

Zheltovodov  et  al  (1990) 

8° 

2.95 

75 

X 

10^ 

2430 

Zheltovodov  et  al  (1990) 

8° 

2.8 

110 

X 

10® 

3950 

All  cases  are  adiabatic  wall 


David  (1993)  performed  the  first  LES  of  a  compression  corner  which  successfully  reproduced  the  Taylor- 
Gortler  vortices  downstream  of  the  shock.  Nevertheless,  the  use  of  a  pseudo-compressible  subgrid-model  did 
not  permit  accurate  quantitative  results.  The  second  and  most  recent  LES  was  Hunt  and  Nixon  (1995)  who 
investigated  the  role  played  by  turbulence,  and  showed  a  direct  correlation  between  the  shock  motion  and 
the  incoming  velocity  fluctuations.  They  also  demonstrated  that  the  size  of  the  separation  bubble  has,  to 
some  extent,  a  weak  effect  on  the  shock  motion.  Despite  the  lack  of  detail  in  the  inner  layer  (a  log-law  wall 
function  was  used  on  a  rough  grid  resolution),  it  displayed  the  qualitative  features  of  the  shock  oscillation 
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observed  experimentally  (Dolling  and  Or  1983). 


A  computation  of  an  adiabatic  turbulent  boundary  layer  flow  past  a  8®  compression  corner  at  Mach  3.0  and 
Re^  =  2  X  10^  was  performed.  As  indicated  in  Table  4,  the  configuration  corresponds  to  the  experiments  of 
Zheltovodov  and  his  colleagues  except  that  the  Reynolds  number  is  lower  for  the  LES.  Table  4  also  includes 
the  flow  conditions  for  other  experiments  and  DNS  which  are  referencec 
LES  results. 

Y 

The  computational  domain  is  -6.0<5  <  x  <  6.05  where  8  is  the  _ 

incoming  boundary  layer  thickness  and  x  is  measured  from  the  , 
corner  (Fig.  24).  The  height  and  width  of  the  computational 
domain  are  Ly  =  3.45  and  Lz  =  4.45  (Le.,  the  same  lengths  as 
previous  flat  plate  configuration).  The  grid  is  1,600,000  tetrahe¬ 
dral  cells.  The  nodal  grid  is  stretched  in  the  y  using  a  geometric  -6-4^0*48 

factor  of  1.088.  The  height  An  of  the  first  cell  adjacent  to  the  tt  i.  i.*  i  j 

1  j  *1  I -i  24  C/0m]pii  vcitrioiisil  domsiixi 

boundary  is  less  than  one  wall  unit.  ®  ^ 

Inflow  data  are  extracted  from  the  previous  flat  plate  boundary  layer  LES  (Fig.  18).  At  each  time  step, 

the  velocity,  pressure  and  temperature  fields  at  the  downstream  station  are  rescaled  and  stored  in  a  datafile. 

Storage  is  performed  over  a  period  equal  to  ZLx/Uoo  (equivalent  to  365/C/’oo)j  which  yields  about  14  GWords. 


Fig.  25  Shape  factor  Fig.  26  Mean  temperature  vs  velocity 
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Fig.  27  Mean  wall  pressure  Fig.  28  Mean  wall  pressure 

The  incompressible  shape  factor  Hi  is  displayed  in  Fig.  25  together  with  the  corresponding  result  for  the 
supersonic  flat  plate  boundary  layer.  The  initial  evolution  is  identical  as  expected,  indicating  that  the 


upstream  boundary  is  located  sufficiently  far  from  the  corner  to  simulate  an  incoming  equilibrium  two 
dimensional  turbulent  boundary  layer  flow. 

In  Fig.  26  the  temperature- velocity  relationship  is  evaluated  at  different  stations  upstream  and  downstream 
of  the  shock.  The  simulation  data  fit  well  with  the  classical  quadratic  profile  originally  due  to  Walz  (1969). 

Fig.  27  displays  the  mean  surface  pressure  on  the  wall.  The  value  of  Re^  for  the  LES  and  the  experiments  of 
Zheltovodov  are  sufficiently  close  that  the  Reynolds  number  effect  should  not  be  significant,  and  therefore, 
since  the  mean  flow  is  unseparated,  scaling  the  distance  x  by  the  boundary  layer  thickness  6  is  reasonable. 
Downstream  of  the  shock,  the  value  of  the  surface  pressure  reaches  the  level  predicted  by  the  experiments 
and  inviscid  theory.  Additionally,  the  slope  of  the  pressure  rise  seems  to  be  correctly  predicted. 


A  indirect  check  can  be  also  be  performed  on  the  wall  pressure,  A  compression  corner  at  the  same  Mach 
number  and  Reynolds  number  but  higher  ramp  angle  a  should  provide  a  smoother  pressure  gradient,  as 
observed  experimentally  (Zheltovodov  ei  al  1990).  The  compression  corner  DNS  from  Adams  (1997)  uses 
the  same  Mach  number,  nearly  the  same  Reynolds  number,  and  a  different  corner  angle  a  =  18®  (Table 
4).  Fig.  28  compares  surface  pressure  along  the  wall.  The  data  are  here  normalized  with  the  downstream 
pressure  level  (maximum  pressure)  in  order  to  take  into  account  the  difference  in  shock  strength.  The  LES 
using  a  =  8®  effectively  predicts  a  sharper  slope  than  the  DNS  using  a  =  18®. 


Fig.  29  Mean  velocity 
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Fig.  31  Reynolds  normal  stress 
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Fig.  30  Mean  temperature 
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Fig.  32  Reynolds  shear  stress 


The  mean  velocity  eX  x/6  =  —3.3  and  3.5  is  shown  in  Fig.  29.  The  computed  profile  displays  close  agreement 
with  the  experimental  data  of  Zheltovodov  et  aL  The  decrecise  in  velocity  due  to  the  shock  is  well  reproduced 
in  the  computation  except  in  the  immediate  vicinity  of  the  wall  where  the  difference  between  the  computed 
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and  experimental  Reynolds  number  is  likely  to  have  an  effect.  The  mean  static  temperature  at  xjb  =  —3.3 
and  3.5  is  shown  in  Fig.  30.  The  agreement  between  computation  and  experiment  is  very  good. 

The  mean  streamwise  resolved  turbulent  kinematic  normal  stress  ^  normalized  using  the  local 

mean  density  *C  ^  and  upstream  wall  shear  stress  r^e/,  is  shown  in  Fig.  31.  The  reference  wall  shear 
stress  is  constant  and  measured  at  x  =  —26,  Comparisons  are  presently  shown  with  Muck  et  al  (1984) 
(a  =  8^;  M  =  3;  Re^  =  1.6  x  10^).  In  Fig.  32  the  Reynolds  shear  stress  distributions  are  shown  for  the  same 
experiment.  The  deformation  of  these  quantities  associated  with  the  unfavorable  pressure  gradient  is  evident 
in  the  vicinity  of  the  corner.  Their  evolutions  show  a  qualitative  agreement.  The  Reynolds  normal  stress  is 
amplified  by  a  factor  as  high  as  4.  This  maximum  increase  appears  at  y/6  =  0.3,  which  is  consistent  with 
experiment.  The  similarity  of  the  experimental  data  between  0.6  <  y/6  <  1.0  for  different  x  is  also  evident  in 
the  LES.  Similar  tendencies  and  qualitative  agreement  with  the  computations  are  found  with  the  Reynolds 
normal  stress  data  of  Zheltovodov  and  Yakovlev  (1984,  1986),  and  Zheltovodov  et  al  (1990).  It  is  likely 
linked  with  the  deformation  of  the  large  scale  structures  in  the  boundary  layer.  No  conclusion  can  be  drawn 
for  the  evolution  of  the  typical  peak  in  the  near  region  due  to  its  lack  even  in  upstream  experimental  data. 
Essentially  the  same  holds  for  the  Reynolds  shear  stresses.  Their  amplifications  correspond  to  experimental 
predictions. 

rau  T 


Fig.  33  Instantaneous  density  at  z  =  2.26  Fig.  34  Instantaneous  temperature  at  z  =  2,26 


Figs.  33  and  34  show  the  instantaneous  density  and  temperature  contours  within  the  boundary  layer  in  an 
X  —  y  plane.  Large  turbulent  structures  can  clearly  be  seen,  both  upstream  and  downstream  of  the  corner 
where  the  shock  appears. 
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